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Abstract. This is the second in a series of papers on the numerical treatment 
of hyperelliptic theta-functions with spectral methods. A code for the numer- 
ical evaluation of solutions to the Ernst equation on hyperelliptic surfaces 
of genus 2 is extended to arbitrary genus and general position of the branch 
points. The use of spectral approximations allows for an efficient calculation of 
all characteristic quantities of the Riemann surface with high precision even in 
almost degenerate situations as in the solitonic limit where the branch points 
coincide pairwise. As an example we consider hyperelliptic solutions to the 
Kadomtsev-Petviashvili and the Korteweg-de Vries equation. Tests of the nu- 
merics using identities for periods on the Riemann surface and the differential 
equations are performed. It is shown that an accuracy of the order of machine 
precision can be achieved. 



1. Introduction 

Solutions to integrable differential equations in terms of theta-functions were intro- 
duced at the beginning of the seventies, see [I. for an account of the history. In 
general, they describe periodic or quasi-periodic solutions. In contrast to the well- 
known solitonic solutions in terms of elementary functions, the theta-functions are 
special transcendental functions defined on some Riemann surface. All quantities 
entering the solution are given in explicit form via integrals which depend implicitly 
on the branch points of the Riemann surface. A full understanding of the functional 
dependence on these parameters seems to be only possible numerically. Algorithms 
have been developed to establish such relations for rather general Riemann surfaces 
as in 1201 or via Schottky uniformization (see QIE]), which have been incorporated 
successively in numerical and symbolic codes, see 0] 00 El E3 an d references 
therein (the first two references are distributed along with Maple 6, respectively 
Maple 8, and as a Java implementation |21|). 

These codes are convenient to study theta-functional solutions on a given surface. 
To establish the dependence of the characteristic quantities of a Riemann surface 
on the branch points, an efficient calculation of the periods is necessary. This is 
especially true if the modular dependence of the theta-functions is of interest as in 
the case of theta-functional solutions to the Ernst equation [Hj in |15l ITT)| . 

In ^21) henceforth referred to as I, we have presented a code for the numerical 
treatment of hyperelliptic solutions to the Ernst equation on a surface of genus 2. 
The integrals entering the solution are calculated by expanding the integrands as a 
series of Chebyshev polynomials using a Fast Cosine Transformation in MATLAB 
and then integrating the polynomials in an appropriate way. In the present arti- 
cle, this code is extended to hyperelliptic surfaces of arbitrary genus and of branch 
points in general position. As an example we consider hyperelliptic solutions to the 
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Kadomtsev-Petviashvili (KP) equation and the Korteweg-de Vries (KdV) equation. 
The precision of the numerical evaluation is tested by checking identities for peri- 
ods on Riemann surfaces and the differential equations. We show that an accuracy 
of the order of machine precision (~ KP 14 ) can be achieved for branch points in 
general position with 32 polynomials and in the case of almost degenerate surfaces 
which occurs e.g. in the well-known solitonic limit with at most 256 polynomials. 
Consequently the solitonic limit can be carried out numerically with machine pre- 
cision which will be shown at the example of the 2-soliton solution to the KdV 
equation. This makes it possible to study the parameter space for hyperelliptic 
surfaces numerically. 

This paper is a sequel to I which was devoted to genus 2 solutions to the Ernst 
equation. We will only briefly repeat the methods already outlined there and refer 
the reader to I for details. The article is organized as follows: in section we 
collect useful facts on the KP and the KdV equation and hyperelliptic Riemann 
surfaces, in section [3] we summarize basic features of spectral methods and explain 
our implementation of various quantities. The calculation of the periods of the 
hyperelliptic surface is performed together with tests of the precision of the nu- 
merics. The theta-series is approximated as in [5] as a finite sum. It is checked 
numerically how well the theta-functional solution solves the differential equation. 
In section 0] we present several examples of hyperelliptic solutions to the KP and 
the KdV equation in general cases and in almost solitonic situations. In section [5] 
we add some concluding remarks. 

2. KP and KdV equation and hyperelliptic Riemann surfaces 

The KP equation for the real valued potential u depending on the three real coor- 
dinates (cc, y, t) can be written in the form 

(1) 3u yy + d x (6uu x + u xxx - 4ut) = 0. 

The completely integrable equation has a physical interpretation as describing the 
propagation of weakly two-dimensional waves of small amplitude in shallow water 
as well as similar physical processes, see [Q. We note that there is a variant of the 
KP equation having a different sign of the u xxx -term which will not be considered 
here. 

Algebro-geometric solutions to the KP equation can be given on an arbitrary Rie- 
mann surface. Here we will only consider hyperelliptic surfaces S g of genus g 
corresponding to the plane algebraic curve 



where we have defined the symmetric (in the branch points) functions S n , n = 
1, . . . , 2g + 2. We will concentrate on real surfaces with real branch points Ei, Fi. 
The branch points are ordered, E\ < F\ < ... < E g+ i < F g+ i. 

Hyperelliptic Riemann surfaces are important since they show up in the context 
of algebro-geometric solutions of various integrable equations such as KdV, Sine- 
Gordon and Ernst. Whereas it is a non-trivial problem to find a basis for the 
holomorphic differentials on general surfaces (see e.g. it is given in the hyper- 
elliptic case (see e.g. pQ) by 



5+1 



(2) M 2 = l[(K - Ei){K - Fi) =: K 2g+2 - S x K 2g+1 + S 2 K 2s + . . . , 



(3) 
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which is the main simplification in the use of these surfaces. We introduce on £ g 
a canonical basis of cycles (<!&,&&), k = 1, ... , n as in Fig. ^ The holomorphic 
differentials dujk are another basis in the space of holomorphic differentials which 
is normalized by the condition on the a-periods 

(4) / duj k = 2TTi6ik- 

J ai 

This fixes the dujk uniquely given the system of cycles. The corresponding matrix 
of 6-periods is given by Bj/j = J b dujk- The matrix B is a so-called Riemann matrix, 
i.e. it is symmetric and has a negative definite real part. 




Figure 1. Canonical cycles. 



The surfaces we consider here are real, i.e. they have an anti-holomorphic involution 
t : Yj g — > E g with r 2 = 1 . The basis of the homology is chosen to have the property 

(5) T bj = bj, TCtj = -aj, j = l,...,g. 

This implies for the normalized holomorphic differentials du>j the reality condition 

(6) rdojj = dujj. 
The matrix of 6-periods is real in this case. 

The theta-function with characteristics corresponding to the curve £ s is given by 

(7) 6 pq (x|B)= J2 cxp[i(B(p + n),(p + n)) + (p + n,27riq + x)|, 

where x £ C 9 is the argument and p, q € C 9 are the characteristics. We will 
only consider half-integer characteristics in the following. The theta-function with 
characteristics is, up to an exponential factor, equivalent to the theta-function with 
zero characteristic (this is the Riemann theta-function, denoted with 0) and shifted 
argument, 

(8) e pq (x|B) = e(x + Bp + 27riq)exp|i(Bp,p) + (p,27riq + x) 
It has the periodicity properties 

(9) e pq (z + 27rie j ) = e 2 ^e pq (z), 9 pq (z + Be,-) = e - 2m ^-3 B «e pq ( z ), 

where ej is the (/-dimensional vector consisting of zeros except for a 1 in jth position. 

The differentials of the second kind needed here have a pole of order I = 1, 2, 3 at 
infinity. They will be denoted by dQi . The Laurent expansion of the corresponding 
integrals is of the form 

(10) / dVli = k l + o(l), 
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where k is the local parameter in the vicinity of infinity. The singularity structure 
characterizes the differentials only up to an arbitrary linear combination of holo- 
morphic differentials. This arbitrariness is eliminated by imposing the condition 
that all a-periods vanish. The normalized differentials satisfy the reality conditions 
rdtli — dili. The vectors of 6-periods of these differentials are denoted by 



(11) u n = / dn u v n = / rffi 2 , w n = / dn 3 . 

Jb n Jb n Jb n 

The periods are real under the above conditions. For later use we define the constant 
c via the expansion of the integral 0.x, 

f p c 

(12) / rffii = k - - +0(k- 2 ). 

Joo k 

The differentials of the second kind dfli are given explicitly on £„ by 



-l 



(13) cW, = — * dK, 

where the q are determined by the condition of vanishing a-periods, 

(H) dn 2 = 1 ^ 2 dK + ..., 

and 
(15) 

3 - - (IS? - \S 2 )K^ - {-LSI ~ \S X S 2 + \S 3 )Kv 

dQs = 2 W) + ' ' 

where . . . denotes the holomorphic differentials leading to vanishing a-periods. For 
the constant c, we get 

(16) C = T~T + C1 - 

It is well known that the 6-periods of normalized differentials of the second kind can 
be expressed in terms of expansions of the holomorphic differentials in the vicinity 
of the singularity (see e.g. [Q). We write the holomorphic differentials in the vicinity 
of some point a in the form 

00 t n 

(17) duj = } v„— -dt, 

n=0 

where t is the local parameter in the vicinity of a. If we choose a = oo + (the infinite 
point in the upper sheet), we have 

(18) v = U, vi = V, v 2 = W/2. 

Thus, the periods U, V and W can be calculated as the 6-periods of the differentials 
d£li, i = 1,2, 3, or via an expansion of the holomorphic differentials. 

Solutions to the KP equation on the above Ricmann surfaces are given by the 
generalization of the Its-Matveev formula for the KdV equation (see e.g. |17p 

(19) u = 2^1n6(Ua; + Vy + Wt + D)+2c, 

where D e C 9 is an arbitrary real vector. Because of |(HJ and the second logarithmic 
derivative in 1)190. the vector D can always be absorbed in the form of a real 
characteristic. By a coordinate change x — * x + tc/3 one can change a solution u 
to the KP equation by —2c. 
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The reality of the argument of the theta-function implies that the theta-function 
which is an entire function does not vanish. This can be seen from the following 
argument for the elliptic case: 

(20) = 1 + 2 ^2 ex P ( \ m2B j cos H mx ) > 0- 

A generalization of the argument to higher genus is straightforward. 

Probably the most elegant way to show that l|19[) provides a solution to the KP 
equation is the use of Fay's trisecant identity [5]. This is an identity for theta- 
functions which holds for four arbitrary points on a Riemann surface, see |17| for 
details. Here we need the identity in the limit that all four points coincide which 
leads to an identity for derivatives of theta-functions. We introduce the directional 
derivatives on C 9 (V acts on z £ C 9 ) 

(21) £> a = v -V, £>; = Vl -V, DZ = v 2 -V. 
Then for an arbitrary vector z the following identity holds: 

D 4 a In 6(z) + 6(D 2 a In 6(z)) 2 + 3D' a D' a In 0(z) - W a D" a In 6(z) 

-24CiL>2i n e( z ) + 12(10C 2 - 3C?) = 0. 



Here the constants C\, C 2 turn up in the Taylor expansion of the differential d£t\ 
from above (a = oo + ), 

(23) dQi = ~ + 2Ci - (6C 2 - 12C 2 )t 2 + .... 

We put u = 2(D 2 ln9(z) — 2Ci) and get for relation (|22J) after differentiation with 
Dl 

(24) D a {Dlu + 6uD a u - 2D" a u) + 3D' a D' a u = 0. 

Because of (|18|) this equation is for z = XJx + Yy + Wt + D equivalent to the KP 
equation. 

2.1. Reduction to KdV equation and solitonic limit. If the hyperelliptic 
surface E 9 is branched at infinity, i.e. if F g+ \ — * 00 in lf2"jl. the integral ^2 is 
single valued. Thus, all its a- and 6-periods vanish which implies that there is 
no y-dependence in the formula (flfjf) . Then the KP equation reduces to the KdV 
equation, 

(25) 6uu x + u xxx - Aut = 0. 

The hyperelliptic surface is now defined by the algebraic curve 

9 

(26) = (K-E g+1 ) Y[(K-Ei)(K-F t ) =: K 2 3+ 1 -S 1 K 29 + S 2 K 2 °-. . .+S 2g+1 . 

1=1 

All definitions apply as before. We use a cut-system of the form Fig. ^ where all 
6-cycles start at the cut [E g+ \ 1 00] along the positive real axis. 

As a local parameter at infinity we can use iVA. The differentials of the second 
kind we need here read 

(27) d0l = _^_ (AS + ClA9 -i + ... Cs)) 
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and 

3ldX ( X g+1 _ I 

2^(A) 1 2 



(28) dn 3 = -^- 7 ^(\9+ 1 - -S^ + d^- 1 + ... + d g ), 



where the constants Cj and di are again defined by the condition that all a-periods 
of the above differentials vanish. Writing the normalized holomorphic differentials 
as 

9 \g-k 

(29) du n = Y dX, 



we get for the 6-periods of the differentials of the second kind 

Si 
2 



Si 

(30) U n = c„i, W n =c n2 + —c nl . 



The constant c reads 

(31) C = -(y +Ci 

An interesting limiting case of the theta-functional solutions on a genus g surface 
is the so-called solitonic limit, see E|- I n this case Ei — ► Fj for i = 1, . . . , g 
which leads to the g-soliton solution. Since g of the cuts collapse to double points, 
the diagonal elements of the Riemann matrix diverge in the used cut-system as 
Bn ~ 21n(.Fj — Ei) whereas all other a- and 6-periods remain finite. The theta- 
series thus breaks down to a sum containing only elementary functions. To obtain 
the standard form of the g-soliton, we choose the vector D in l|19|) to correspond 
to the half-integer characteristic 



(32) \ 



\ ... \ 
... 



On an elliptic surface we get with the above relations for (|f 9[) the 1-soliton solution 
of KdV equation, 

U 2 

(33) u = J— + 2c ' z = Ux + wt > 

2 cosh | 



where 



(34) U = 2 y /E 2 -E 1 , W = -2^E 2 - E X S U c = 

Similarly we get for the 2-soliton 

U\ - Ul Ul cosh 2 f + E/f sinh 2 f 



2 ' 



(35) u - 2c = 



where 



(£/i cosh §■ cosh f-U 2 sinh sinh f) 2 



(36) Zl = U iX + Wit, Ut = 2 y /E 3 -E 1 , U 2 = 2^E 3 

and 



(37) Wi = -^E 3 -E 1 {E 3 + 2EJ, W 2 = -^E 3 -E 2 {E 3 + 2E 2 



E 3 
' 2 ' 



A generalization to higher genus is straight forward, see ^H]. Formulas H34|) and 
(|35|l hold also for the KP equation with the appropriate definition of z and U. 
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3. Numerical implementations 

The numerical task in this work is to approximate and evaluate analytically defined 
functions as accurately and efficiently as possible. To this end it is advantageous to 
use (pseudo-)spectral methods which are distinguished by their excellent approxi- 
mation properties when applied to smooth functions. Here the functions are known 
to be analytic except for isolated points. In this section we explain the basic ideas 
behind the use of spectral methods and describe in detail how the theta-functions 
and its derivatives can be obtained to a high degree of accuracy. 

3.1. Spectral approximation. The basic idea of spectral methods is to approxi- 
mate a given function / globally on its domain of definition by a linear combination 

N 
fe=0 

where the functions <pk are taken from some class of functions which is chosen 
appropriately for the problem at hand. The spectral coefficients are determined by 
a so-called collocation method, i.e. by the condition that f(xi) = XifeLo a k4>k(xi) at 
selected points (xi)i = o : n- Since we are interested in an approximation of functions 
defined on a finite interval, orthogonal polynomials, in particular Chebyshev T n {x) 
polynomials will be chosen. They are defined on the interval I = [—1,1] by the 
relation 

T„(cos(f)) = cos(nt), where x = cos(t), t € [0,7r]. 
The Chebyshev polynomials are used because they have very good approximation 
properties and because one can use fast transform methods when computing the 
expansion coefficients from the function values provided one chooses the collocation 
points xi = cos(irl/N) (see ^U] and references therein). 

Let us briefly summarize here some basic properties of Chebyshev polynomials: 
The addition theorems for sine and cosine imply the recursion relations 

(38) ~~n~+Y~ = 2Tn{x) 

for their derivatives. The Chebyshev polynomials are orthogonal on / with respect 
to the hermitian inner product 

(f,9) = J J{x)g{x)-j==. 

We have 

(39) (T m ,T n ) c m — S mn 
where cq = 2 and c; = 1 otherwise. 

The fact that / is approximated globally by a finite sum of polynomials allows us to 
express any operation applied to / approximately in terms of the coefficients. Let us 
illustrate this in the case of integration. So we assume that / = pn — Y^,n=o a nT n , 
and we want to find an approximation of the integral for pjy, i.e., the function 



F(x)= / f(s)ds, 



l 

so that F'(x) — f(x). We make the ansatz F(x) = ^2n=obnT n (x) and obtain the 
equation 

JV N 
n=0 n=0 



F' = Vfc„r; = Va„r„ = /. 
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Expressing T n in terms of the T' n using l|38[) and comparing coefficients determines 
all 6; in terms of the a n except for bo . This free constant is determined by the 
requirement that F(—l) = 0. Thus, to find an approximation of the definite integral 
J_j f(x)dx we proceed as described above, first computing the coefficients a n of /, 
computing the b n and then calculating the sum of the necessary coefficients. 

3.2. Numerical treatment of the periods. The quantities entering formula Ijl9(l 
are the periods of certain differentials on the Riemann surface. The value of the 
theta-function is then approximated by a finite sum. 

The periods of a hyperelliptic Riemann surface can be expressed as integrals be- 
tween branch points. Since we need in our example the periods of the holomorphic 
differentials and differentials of the second kind with poles at oo , we have to 
consider integrals of the form 

f p i K n dK 

(40) I — , n = 0,..., 5 + 3, 

where the Pj, j = 1, . . . , 2g + 2 denote the branch points of 

We parametrize the straight line between the two branch points Pi and Pj by 

K = \{P 3 +P l ) + \(P J -Pd 
to obtain the integral H40(l in the normal form 

(41) f Q0 + Qlf+ /-- + Q ^ 9+3 Hit) dt, 
J-i VI — t 2 

where the a, are complex constants and where H(t) is a continuous (in fact, an- 
alytic) complex valued function on the interval [—1,1]. This form of the integral 
suggests to express the powers t n in the numerator in terms of Chebyshev polyno- 
mials and to approximate the function H(t) by a linear combination of Chebyshev 
polynomials 

n>0 

The integral is then calculated with the help of the orthogonality relation Q39J1 of 
the Chebyshev polynomials. 

This is the procedure used in I for general position of the branch points which can 
also be applied here for low genus. The problem is that there does not seem to 
be a direct way to express the terms K n in H40(l after the linear transformation in 
terms of Chebyshev polynomials in t. This is an obstacle to the determination of 
the periods for general genus. Another problem are situations close to the solitonic 
limit. Since the cut-system Fig. ^ is adapted to this case, the a-periods do not lead 
to problems. For the 6-periods, the function H in the expression for bi B41JI behaves 
like 1/y/t + Ei — Fi near t = 0. For Ei ~ Fi this behavior is only satisfactorily 
approximated by a large number of polynomials. The use of a large number of 
polynomials however does not only require more computational resources but has 
the additional disadvantage of introducing larger rounding errors. It is therefore 
worthwhile to reformulate the problem since a high number of polynomials would 
be necessary to obtain optimal accuracy in this case. 

It is possible to address both problems in one approach. The idea is to use sub- 
stitutions in the integrals (|40|l leading to a regular integrand. To determine the 
a-periods, we use 

(42) K = ^--^ + £l— : =i cosh x. 
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After a linear transformation which transforms the integration path to the interval 
[— 1, 1], the integral is computed with the Chebyshev integration routine sketched 
above. This also works in situations close to the solitonic limit. To treat the b- 
periods in this case, we split the integral from Fi to Ei + \ in two integrals from Fi 
to (Fi + E i+ i)/2 and from (Fi + E i+ i)/2 to E i+ i. In the former case we use the 
substitution (|42fl and 

(43) E i+1 +F i+1 + Fi+i-E i+ i 

in the latter. These substitutions lead to a regular integrand even in situations close 
to the solitonic limit. After a linear transformation, the integrals are computed with 
the Chebyshev integration routine. The disadvantage of this approach is that the 
definition of the square root adapted to the cut-system in I via monomials cannot 
be used in the case of the above non-linear transformations. Here the root is chosen 
in a way that the integrand is a continuous function on the path of integration. This 
can be achieved in general for real branch points. 

To test the numerics we use the fact that in our case the integral of any holomorphic 
differential along a contour surrounding the cut [-Ei,-Fi] in positive direction is 
equal to minus the sum of all a-periods of this integral. Since this condition is not 
implemented in the code it provides a strong test for the numerics. The symmetry 
of the Riemann matrix (the matrix of o-periods of the holomorphic differentials 
after normalization) can be used to estimate the error in the numerical evaluation 
of the o-periods. We define the function err as the maximum of the norm of the 
difference in the a-periods discussed above and the difference of the off-diagonal 
elements of the Riemann matrix. Additional tests for the periods U, V and W 
follow from (|18Jl . We show in Fig. ||a genus 2 situation with the branch points 
[— 3, — 2, — 1, — 1 + + e]. It can be seen that 32 polynomials are sufficient 

in general position (e — 1) to achieve optimal accuracy. Since MATLAB works 
with 16 digits, machine precision is in general limited to 14 digits due to rounding 
errors. These rounding errors are also the reason why the accuracy drops slightly 
when a higher number of polynomials is used. In an almost degenerate situation 
(e = 10~ 14 ), optimal accuracy is achieved with 128 polynomials. Note that all 
quantities entering (|19fl are within 10 -14 and better equal to the limits given in 
(|36|l and (|37|l . Thus the solitonic limit can be reached numerically with machine 
precision. 

The above procedure to determine the periods can be directly implemented within 
MATLAB for arbitrary genus. Due to the efficient vectorization algorithms of 
MATLAB, the calculation of the periods with optimal accuracy takes only a sec- 
ond even for genus 10 on the used low-end computers. The limiting factor is here 
whether the matrix A of a-periods is ill conditioned. This is the case if a consid- 
erable number of the entries of A are of the order of the rounding error (10 -16 ). 
Thus due to the limited number of digits, the inversion of the matrix A which is 
necessary to determine the Riemann matrix can only be carried out with reduced 
accuracy which is independent of the number of polynomials used in the spectral 
expansion. For large genus, this becomes the limiting factor in the determination of 
the Riemann matrix. For instance in the case of the genus 20 surface with branch 
points [—21, 20, . . . , 20], the identity for a-periods is satisfied with 32 polynomials 
up to 10 -15 , the symmetry of the Riemann matrix only up to 10 -6 since the matrix 
A is badly scaled. The calculation of the periods takes less than a second in this 
case. 



10 



J. FRAUENDIENER AND C. KLEIN 



A e=1 
□ e=10~ 



6 7 
log N 



Figure 2 . Test of the numerics for the a-periods at several points 
in the space-time. The error is shown in dependence of the number 
N of Chebyshev polynomials. 



3.3. Theta-functions. The calculation of the theta-functions below is applicable 
for arbitrary Riemann surfaces, it is not limited to the hyperelliptic case. The 
theta-series (J7J) for the Riemann theta-function (the theta-function in (0) with zero 
characteristic, theta-functions with characteristic follow from (JHJ) is approximated 
as the sum 

N N f 

(44) 6(x|B)= E ■■■ E expj-(M,BM) + (M,x) 

where M is a vector in Z 9 with the components m, . . . , n g . We use the periodicity 
properties of the theta-function <0 to minimize the absolute value of x. The value 
of N is determined by the condition that terms in the series (0 for n > N are 
strictly smaller than some threshold value e which is taken to be of the order of 
1CP 16 . To this end we determine the eigenvalues of B and demand that 

(45) N > (N| + VINI 2 + 2]negB max ) , 

where B max is the real part of the eigenvalue with maximal real part (B is nega- 
tive definite). Similar formulas can be obtained for theta-derivatives. For a more 
sophisticated analysis of theta summations see |S] . In the studied examples we 
found values of N between 2 and 40. If the eigenvalues of B differ by more than an 
order of magnitude which can happen close to partial degeneration of the surface, a 
summation over an ellipse rather than over a sphere in the hypercube (rii, . . . , n^g), 
i.e. different limiting values for the rij as in 5 could be more efficient. The summa- 
tion over the hypercube has, however, the advantage that it can be implemented 
in MATLAB for arbitrary genus. In addition it makes full use of MATLAB's vec- 
torization algorithms outlined below. Thus it is questionable whether a summation 
over an ellipse would be more efficient in terms of computation time in this setting. 

Here we made use of MATLAB's efficient way to handle matrices. We generate a 
2N + 1-dimensional array containing all possible index combinations and thus all 
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components in the sum l|44ll which is then summed. To illustrate this we consider 
the simple example of genus 2 with N = 2. The summation indices are written as 
(2N + 1) x (27V + l)-matrices since g = 2. Each of these matrices contains 2N + 1 
copies of the vector with integers — (2N + 1), . . . , 2N + 1. N 2 is the transposed 
matrix of N\. Explicitly, we have 



(46) Ni 



The terms in the sum Q44H can thus be written in matrix form 



/ 


2 


2 


2 


2 


2 


\ 




( 2 


1 





-1 


-2 \ 




1 


1 


1 


1 


1 






2 


1 





-1 


-2 





















, N 2 = 


2 


1 





-1 


-2 




-1 


-1 


-1 


-1 


-1 






2 


1 





-1 


-2 


V 


-2 


-2 


-2 


-2 


-2 


/ 




V 2 


1 





-1 


- 2 / 



(47) exp(~Ni*JVi.Bn 



Ni * N 2 B 12 



-N 2 * N 2 B 22 + Nx * xi + N 2 * x 2 



where the operation N\-kN 2 denotes that each element of N\ is multiplied with the 
corresponding element of N 2 . Thus, the argument of exp is a (2N + 1) x (2N + 1)- 
dimcnsional matrix. Furthermore, the exponential function is understood to act 
not on the matrix but on each of its elements individually, producing a matrix of 
the same size. The approximate value of the theta-function is then obtained by 
summing up all the elements in (|47|l . 

The most time consuming operations are the determination of the bilinear terms 
involving the Riemann matrix. If one wants to calculate solutions to the KP equa- 
tion, these terms only have to be determined once for a given surface. The integer 
N is in this case fixed for the largest ||x|| in the plot. We note that the summation 
is very fast even though the used determination of N is rather crude. For instance 
in the case of genus 2 with N = 100, the calculation of the theta-function takes 
0.1s on the used computers. 

The limiting factor is here the available memory since arrays of the order (27V + l) 9 
have to be multiplied with each other. On the used low-end computers we could 
deal with rather general genus 4 situations, but the limit was reached for genus 6 
and N = 5. The summation is still very efficient, the calculation of the bilinear 
terms and the determination of the coefficients took 16s in the latter case, the 
subsequent calculation of the linear terms in l|47|l and the summation, which have 
to be carried out for each value of x and t, took roughly 4s. Thus the limitations 
we had to face were not due to the computing time but due to missing memory. 

The quality of the theta summation can be checked directly by putting the ap- 
proximate solution into the differential equation. The necessary derivatives of the 
theta-function entering the function u and its derivatives can be calculated directly, 
e.g. 

N N f . 

(48) ^6(z,B) = •■■ E ((M,U)rexp|-(M,BM) + (M,z)|, 



m=— N 



n„ = -N 



where TV is determined as TV for the theta-series. The differential equations could 
be satisfied to the order of 10 -13 and better. For the plots in the next section we 
use an accuracy of the order of 10~ 10 . The error is mainly due to the fact that we 
did not choose N big enough that also the terms in 6th derivative of the form l|48|l 
are of the order of the rounding error. For the solution only the second derivative 
is important. The accuracy of the solution is thus in general much higher as can be 
seen for instance in the solitonic example where the differential equation is satisfied 
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to the order of 10 13 , but the difference to the solitonic solution is smaller than 

io- 15 . 



In this section we present plots of hyperelliptic solutions to the KdV and the KP 
equation as examples for the previously discussed code. Such plots can already 
be found in ^| |2] and We will show general situations as well as almost 
degenerate surfaces which are identical to their corresponding solitonic solution up 
to numerical accuracy. In all plots, the vector D in l|19f> is chosen to correspond to 
the characteristic Ij32l) . 

It turns out that almost solitonic solutions take less computational resources. The 
calculation of the periods requires the use of more polynomials, but the most time- 
consuming part is the theta summation. Since the diagonal elements of the Riemann 
matrix diverge in this limit, less terms have to be considered in the summation. 
The generation of the plots is thus considerably faster than in general position of 
the branch points. 

To begin we want to discuss plots of genus 2 solutions to the KdV equation. We 
consider a hyperelliptic surface of the form l|26l) with branch points [—2,-2 + 
e, —1, —1 + e, 3]. In Fig. [3] we show the case for e = 10~ 14 which is identical to the 
2-soliton solution within machine precision. The difference between the shown plot 
and the solution (|35|1 is actually less than 1CV 15 . The one dimensional waves in 
shallow water are depicted in dependence on x and t. It can be seen that a soliton 
coming from the right and traveling in negative x-direction has a collision with a 
soliton traveling in positive ^-direction. At the collision the typical phase shift can 
be observed, otherwise the solitons are unaffected. 



4. Hyperelliptic solutions to KdV and KP equations 
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Figure 3. 'Almost' solitonic genus 2 solution to the KdV equation. 



In Fig. Q]we show the case e = 1. The almost periodic nature of the solution is 
clearly recognizable. The solution can be interpreted as an infinite train of solitons. 
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Figure 4. Almost periodic genus 2 solution to the KdV equation. 

To obtain solutions on a surface of genus 6 we considered the surface with the branch 
points [-6, -6 + e, -4, -4 + e, -2, -2 + e, 0, ep, 2, 2 + e, 4, 4 + e, 6] . The situation for 
e = 1 is shown in Fig.|SJ The depicted situation is at the limit imposed by memory 
on the used computer. In Fig. we show a solution of genus 6 in an almost solitonic 




Figure 5. Almost periodic genus 6 solution to the KdV equation. 

situation. One can see the collision of 6 solitons at the center of the plot. 

A genus 4 solution of the KP equation is shown for fixed t on the hyperelliptic 
surface with branch points [—5, —4, — 3, — 3 + e, — 1, — 1 + e, 1, 1 + e, 3, 3 + e]. In Fig. [7] 
we show the almost periodic situation for e = 1. The almost solitonic case with 
e = 0.001 is shown in Fig.HJ 
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~20 . 



Figure 6. Almost solitonic genus 6 solution to the KdV equation. 




Figure 7. Almost periodic genus 4 solution to the KP equation. 



5. Conclusion 

In this article we have presented a scheme based on spectral methods to treat 
hyperelliptic theta-functions of in principle arbitrary genus numerically. It was 
shown that an accuracy of the order of machine precision could be obtained with an 
efficient code. As shown, spectral methods are very convenient if analytic functions 
are approximated. Close to singularities such as the degeneration of the Riemann 
surface, analytic techniques must be used to obtain analytic integrands as in the 
discussed example. This made it possible to study the solitonic limit numerically 
within machine precision. 
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Figure 8. Almost solitonic genus 4 solution to the KP equation. 

To consider the case that more than 2 branch points coincide, a higher number of 
polynomials or additional analytical techniques have to be used. As shown in I this 
makes it possible to study the parameter space of hyperelliptic surfaces including 
almost degenerate cases with at most 512 polynomials with maximal accuracy (here 
we used at most 128 polynomials). Both the calculation of the periods and the 
theta summation are very efficient. The limiting factor for the periods is the finite 
number of digits (16 in MATLAB) which leads to a drop of accuracy if the matrix 
of a-periods is ill-conditioned. The theta summation is limited by the available 
memory 

The presented numerical code is able to treat general hyperelliptic surfaces. The 
used method can in principle be extended to more general Riemann surfaces if the 
differentials are known explicitly as e.g. in the case of Z^-curves. 
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